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1.  Introduction 


In  previous  work  (Powell  et  al.,  1993;  Powell  and  Zielinski,  1995,  1997,  1999,  2001,  2003, 

2004),  we  have  developed  a  number  of  two-dimensional  (2-D)  models  for  solving  coupled 
Maxwell  and  energy  transport  equations  in  two  dimensions.  The  models  have  been  used  to 
investigate  current  and  heat  transport  in  the  rail  and  annature  of  both  solid-  and  plasma-armature 
railguns.  All  of  these  calculations  have  been  based  upon  the  assumption  that  the  rail  height  h  is 
very  large  when  compared  to  the  rail  separation  w  (see  figure  1).  The  models  are  known  to 
predict  values  of  the  in-bore  fields  and  annature  acceleration  that  are  too  large  when  compared  to 
realistic  experimental  values  for  finite-rail-height  launchers.  Some  methods  for  scaling  the 
results  to  make  them  more  realistic  have  also  been  studied  (Powell  and  Zielinski  1995,  1997, 
2004). 


When  the  rail  height  is  infinitely  large,  there  exists  only  one  component  of  the  magnetic 
induction  (along  the  direction  of  the  rail  height)  and  two  components  of  the  electric  field  and 
current  density  (along  directions  normal  to  the  rail  height).  In  this  report,  we  consider  the 
“alternate”  2-D  model  in  which  the  current  is  conducted  only  in  one  direction,  while  the 
magnetic  induction  has  two  components.  The  model  should  provide  us  with  a  methodology  for 
considering  a  number  of  different  problems  that  were  not  possible  with  our  previous  model. 


1 


Such  problems  include,  for  example,  electromagnetic  (EM)  field  transport  in  finite-height  rails  at 
distances  far  removed  from  the  armature.  At  these  distances,  we  expect  that  current  is  conducted 
only  along  the  direction  of  the  rail,  but  that  the  magnetic  induction  field  will  “wrap  around”  the 
rail.  Consequently,  there  will  exist  only  one  component  of  the  current  density  (along  the  rail 
direction)  but  two  components  of  the  magnetic  induction  (along  directions  normal  to  the  rail 
direction).  We  should  also  point  out  that  we  have  recently  extended  a  similar  model  in 
cylindrical  coordinates  that  was  used  to  investigate  the  disruption  of  shaped  charge  jets  by 
axially  directed  currents.  That  model  has  now  been  used  to  study  current  diffusion  in  some  fairly 
complicated  cylindrical  geometric  configurations  in  another  report  (Powell,  2008). 

In  our  previous  work,  we  solved  Maxwell’s  equations  directly  for  the  EM  fields,  but  in  this 
report  we  employ  the  formalism  of  the  vector  and  scalar  potentials.  This  approach  is  used 
because  there  is  a  single  component  of  the  vector  potential  from  which  the  EM  fields  follow 
easily.  The  calculations  can  therefore  be  significantly  simplified;  the  two  approaches  yield 
similar  results. 

Very  complex  models  and  computer  codes  for  doing  similar  types  of  calculations  have  been 
developed  in  the  past  by  others.  The  most  noteworthy  example  for  application  to  railguns  is 
perhaps  the  simulation  code  EMAP3D  (Electromechanical  Analysis  Program  in  Three 
Dimensions)  (Hsieh,  1995).  That  code,  which  is  three  dimensional,  has  the  capability  of 
undertaking  the  sorts  of  problems  of  interest  to  us.  Our  intent,  however,  is  to  develop  a  simpler 
model  with  which  we  can  undertake  elementary  calculations  and  that  can  be  extended  to  more 
complicated  studies  in  the  future.  Two-dimensional  models  that  are  similar  to  the  one  developed 
here  include  models  by  Kerrisk  (1982),  who  investigated  current  diffusion  in  the  rails  of  a 
railgun,  and  by  Freeman  (1988),  who  studied  the  launch  of  projectiles  with  a  reconnection  gun. 

The  primary  purpose  of  the  report  is  to  develop  the  model  and  to  provide  the  methodology 
necessary  for  undertaking  future  calculations.  One  simple  example,  however,  is  calculated.  The 
example  is  concerned  with  the  diffusion  of  current  and  the  resulting  ohmic  heating  of  two  long 
parallel  rails  that  carry  current  only  along  the  direction  of  the  rails.  This  problem  has  potential 
importance  for  investigating  rail  erosion  and  assessing  effects  that  might  contribute  to  transition. 
We  also  undertake  a  calculation  of  the  inductance  per  unit  length  L'  of  the  pair  of  rails. 

The  organization  of  the  report  is  as  follows.  In  section  2,  we  derive  the  general  governing 
equations  in  three  dimensions  and  then  obtain  the  limiting  results  for  a  2-D,  Cartesian  model. 
Some  of  this  analysis  was  contained  in  our  previous  report  (Powell,  2008),  but  it  is  also  included 
here  for  completeness.  In  section  3,  we  undertake  the  current  diffusion  calculation  just 
described.  Finally,  section  4  contains  a  summary  and  indicates  possible  future  directions  for  this 
research. 
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2.  Model  and  Formalism 


We  consider  one  or  more  electrical  conductors  contained  within  some  domain  D.  The 
conductors  carry  current,  which  may  be  specified  a  priori  or  may  need  to  be  calculated.  The 
conductors  may  be  either  in  a  solid  or  liquid  state,  and  they  may  be  moving  with  a  velocity  v. 
We  do  neglect  internal  motion,  however,  so  that  the  conductors  move  as  perfectly  rigid  bodies. 
Our  intent  is  to  determine  the  EM  fields  throughout  the  domain  as  well  as  the  temperature  of  the 
conductors  as  a  function  of  both  space  and  time. 

2,1  Vector  Potential  Formulation  in  Three  Dimensions 

We  begin  with  Maxwell’s  equations  and  Ohm’s  law  in  three  dimensions.  Let  E,  B  ,  v  ,  and  J 
be  the  electric  field  intensity,  the  magnetic  induction  field,  the  velocity  of  the  conductor,  and  the 
current  density,  respectively.  Let  //  and  <7  be  the  permeability  and  electrical  conductivity, 
respectively.  Then  we  have 


< 

X 

tqi 

II 

1 

CL)|  g;1 
|  GDI 

(1) 

Vx(B/ju)  =  J, 

(2) 

J  =  a  (E  +  v  x  B), 

(3) 

V»B  =  0. 

(4) 

These  equations  are  written  in  the  low-frequency  approximation,  which  applies  whenever  the 
velocity  v  is  small  when  compared  to  the  light  speed  c o,  and  all  relevant  time  scales  in  the 
problem  are  long  when  compared  to  (//rxr(f)  1 .  Both  conditions  are  easily  met  in  the  electrical 

conduction  problems  that  we  consider.  It  is  also  assumed  that  the  relationship  between  the 
magnetic  field  intensity// ,  which  strictly  should  appear  in  equation  2,  and  the  magnetic 
induction  B  is  given  by  H  =  B  /  //,  where  /i  is  assumed  to  be  constant  within  the  material. 

We  now  define  the  vector  potential  A  by  the  relation 

B  =  VxA.  (5) 

With  this  definition,  it  is  apparent  that  equation  4  is  satisfied  automatically,  since  the  divergence 
of  the  curl  of  any  vector  is  always  zero.  Substitution  of  equation  5  into  equation  1  then  provides 
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Vx 


dt  j 


0. 


(6) 


This  equation  has  the  solution 

~E+di =-v*  <7> 

where  Vtp  is  an  arbitrary  function  that  constitutes  the  scalar  potential  contribution  to  the  electric 
field. 

We  finally  observe  that  specification  of  the  curl  of  A  as  given  in  equation  5  does  not  specify  the 
vector  uniquely.  A  unique  specification  can  be  obtained,  however,  if  we  also  indicate  the 
divergence  of  the  vector.  A  convenient  choice  (there  are  others),  which  is  known  as  the 
Coulomb  gauge,  is 

V*A  =  0.  (8) 

If  we  assume  this  condition,  substitute  equations  3,  5,  and  7  into  equation  2,  and  expand,  we 
obtain 

r)T  1  - 

( j  —  =  —  V2A-V(l//t)xVxA  +  <j(v  xV  x  A)  -  oV  (/>,  (9) 

dt  jl 

where  by  definition  V2A  =  V2  A  J  +  V2  A  j  +  V2  A.k.  Equation  9,  solved  with  appropriate 
boundary  conditions  and  a  specification  of  V</>,  permits  determination  of  A  .  Substitution  of 
that  result  in  equations  7,  5,  and  3  then  produces  the  remaining  variables  E,  B,  and  J.  It 
should  also  be  pointed  out  that  equation  8  holds  not  only  in  the  conductors  but  in  the  surrounding 
vacuum  as  well,  provided  we  take  a  =  0. 


2.2  Governing  Equations  for  2-D  Rectangular  Coordinates 

For  the  2-D  case  with  rectangular  coordinates,  we  assume  that  all  variables  depend  only  on  x  and 
y  and  not  on  the  coordinate  z.  Furthermore,  for  the  specific  case  under  study  here,  we  assume 
that  A  =  Ad. ,  E  =  Ed. ,  and  J  =  Ja..  It  is  evident,  then,  from  equation  5  that  v  and  B  may  have 
components  in  both  the  x  and  y  directions.  If  we  furthermore  note  that  equation  8  has  only  a  z 

component,  we  must  have  that  —  =  —  =  0.  Consequently,  the  only  nonzero  component  of 

dx  dv 

can  be  represented  as  E0(t),  where  E0  is  a  position-independent  function.  If  we  now  expand 
equation  9  completely,  we  find  for  the  2-D  equation  governing  the  evolution  of  A  the  result 


IS! 

h 

|S! 

h 

fe 

h 

q 

_  1 

(d2A  d2A 'l 

-L 

i 

/  d/J  dA  djU  dA  ^ 

[dt  dx  dy  j 

M 

l  a*2  a/  j 

v  dx  dx  dy  dy  J 
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In  this  equation,  u  and  w  represent  the  x  and  y  components  of  the  velocity,  respectively.  Once 
equation  10  has  been  solved  for  the  vector  potential  A,  we  can  determine  the  magnetic  induction 
and  current  density  from  equations  5,  3,  and  7.  For  the  2-D  case,  we  have 


and 


D  -dA 

^ x  ’ 
qy 

(ID 

O  _  dA 
y  dx’ 

(12) 

'  dA  dA  dA  ' 

—a - 1-  h,Q  +  it - 1-  w — 

v  dt  dx  dv  j 

(13) 

In  general,  the  electrical  conductivity  cr  depends  sensitively  on  the  temperature  of  the  conductors. 
Consequently,  it  is  necessary  to  solve  an  energy  transport  equation  within  the  conducting  regions 
in  order  to  detennine  the  temperature  as  a  function  of  space  and  time.  We  account  for  ohmic 
heating  and  ordinary  heat  conduction,  and  we  allow  for  melting  but  no  other  phase  change.  Let 
e,  C,  p,  and  k  be  the  specific  internal  energy,  the  specific  heat  capacity,  the  density,  and  the 
thennal  conductivity,  respectively.  Let  Tm  be  the  melting  temperature,  and  let  Hf  be  the  latent 
heat  of  fusion.  We  then  assume  that  the  time  evolvement  of  the  energy  is  given  by 


de  de  de  d2T  d2T  die  dT  die  dT  2  , 

dt  dx  dy  dx2  dy2  dx  dx  dy  dy 


(14) 


The  first  four  terms  on  the  right-hand  side  account  for  heat  conduction,  whereas  the  last  tenn 
accounts  for  ohmic  heating. 

The  relationship  between  e  and  T  must  be  determined  by  an  equation  of  state,  and  we  use  the 
simple  representation 

T 

e  =  j  C{T')dT'  +  HfQ(T  -Tm).  (1 5) 

o 


In  this  equation,  0  is  the  Heaviside  function,  which  is  equal  to  zero  when  T  is  less  than  Tm  and  is 
equal  to  unity  otherwise.  The  equation  clearly  predicts  that  melting  occurs  at  a  well-defined 
temperature  Tm.  Since  the  energy  depends  on  ohmic  heating  and  the  electrical  conductivity 
depends  on  temperature,  there  is  an  obvious  coupling  between  equations  10  and  14. 
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2.3  Boundary  Conditions 

Boundary  conditions  that  must  be  applied  at  the  outer  edges  of  the  domain  29  are  determined  by 
physical  considerations  for  each  problem  considered.  There  are,  however,  additional  boundary 
conditions  that  must  hold  at  the  interface  between  the  conductors  and  the  vacuum.  For  the  vector 
potential  A,  these  conditions  can  be  obtained  by  integrating  equations  2  and  4  across  such  an 
interface  of  thickness  S  and  then  taking  the  limit  as  S  — >  0.  Let  n  be  the  unit  normal  to  the 
interface.  Then,  if  we  assume  that  there  are  no  surface  currents,  the  integrations  yield 

n»{B}  =  0  (16) 


and 


n  x{B  /  ju}  =  0, 


(17) 


where  the  braces  are  used  (in  this  section  only)  to  denote  the  change  in  the  quantity  in  question 
as  the  interface  is  crossed.  Thus,  in  accordance  with  standard  results,  we  have  that  the  normal 
component  of  B  and  the  tangential  component  of  H  are  continuous  at  the  boundary.  In  terms 
of  the  vector  potential  A,  equation  16  implies  that  A  is  continuous  at  the  boundary,  whereas 
equation  17  implies  that 


i  a4|  |  f  i  a4| 

M  d  r  |  ‘(//  dy  j 


=  0. 


(18) 


Thus  if  jx  is  discontinuous  at  the  boundary,  so  too  is  the  nonnal  derivative  of  A. 

Since  heat  conduction  is  accounted  for,  there  must  also  be  boundary  conditions  on  the 
temperature  at  the  conductor  surfaces.  The  appropriate  condition  can  also  be  obtained  by 
integrating  equation  14  across  the  boundary.  We  obtain 

/7«{*Vr}  =  0.  (19) 


Therefore,  since  k  is  zero  in  the  vacuum,  equation  19  produces  the  usual  result  that  the  nonnal 
derivative  of  T  must  vanish  at  a  conductor-vacuum  interface. 


2.4  Generalized  Coordinates 

As  in  previous  work,  we  often  find  it  convenient  to  transform  the  governing  equations  to  a  set  of 
generalized  coordinates.  The  purpose  of  this  transformation  is  to  map  a  complicated  physical 
space  into  a  much  simpler  computational  space  in  which  the  calculations  are  performed.  The 
computational  space  is  chosen  to  be  rectangular.  This  procedure  makes  it  possible  to  consider 
irregularly  shaped  boundaries,  for  example,  in  a  rather  straightforward  manner.  Once  the 
computation  has  been  completed,  the  results  are  then  “transformed  back”  to  the  physical  space 
for  presentation. 
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We  now  outline  the  general  procedure  for  perfonning  the  transformation.  Greater  detail  can  be 
found  in  our  previous  work  (Powell  and  Zielinski,  1997)  or  in  various  textbooks  (Anderson  et  ah, 
1984;  Hoffmann,  1989).  We  denote  the  physical  coordinates  by  (x,y,t)  and  the  computational 
coordinates  by  (<;,?/, r),  and  we  write  quite  generally 

€  =  €(x,y,t) 

ij =  n(x,y,t) 

T  =  t.  (20) 

The  partial  derivatives  that  appear  in  equations  10  and  14  must  be  expressed  in  terms  of  the 
computational  coordinates,  and  the  appropriate  relations  can  be  obtained  from  the  chain  rule. 

We  have  for  any  function  F 

Fx=F&+Fr,Vx 

Fy  =  F&  +  FnJ]y 

F,x  =  FtL  +  FJI„  +  £FU  +  n;Fn„  +  'LZjlJPfn 

Fyy  =  F£yy  +  FnJhy  +  %FU  +  <Fn,  + 

F,=F r-  (21) 

We  have  used  here  (this  section  only)  subscripts  that  are  coordinates  to  denote  partial  derivatives 
with  respect  to  those  coordinates.  Consequently,  r/xx  =d2r//dx2,  F::  =  d2F  / dg2 ,  and  so  forth. 

We  now  substitute  from  equation  21  into  equations  10  and  14  to  produce  in  terms  of  the 
computational  variables 

<j\At  +  u{A^x  +  Antjx)  +  w(A^y  +  Antjy )  +  E0\ 

=  (1  ///)[/!=  (£„  +  ^)  +  A  v  (  Jixx  + 

+  a^I  +  g )  +  Am(j]2x  +  ril)  +  2 tfxnx  +  g,n..)A:  ] 

-  (1  /  // 2  )[(ji£x  +  Mnnx  )(A^X  +  Anrjx ) 

+  +  MnVy )(A£y  +  Ajiy)\  ,  (22) 


p[eT  +  u{e£x  +  enrix)  +  w(e^y  +  etJT]v )] 

=  427(4+4+7;(^+/7vv) 

+  &  +  ) +  Tnn  ^  K  ^ 

+  [(*>£  +  KnTlx  )(r^x  +  TnVx) 

+  (k tig y  +  Knriy)(T^y  +  Tri7iy)  +  J2  !a\  (23) 
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As  is  generally  the  case  in  these  types  of  transformations,  the  transformed  equations  are  much 
more  complicated  than  the  original  ones. 

Our  earlier  work  (Powell  and  Zielinski,  1997,  2001)  contained  both  complicated  and  simple 
transformations  as  indicated  in  equation  20  in  order  to  achieve  some  desired  mapping.  For  the 
problem  to  be  considered  here,  only  the  trivial  transfonnation 

4  =  y 

n  =  x  (24) 

is  required.  Consequently,  we  do  not  discuss  the  details  of  more  complicated  transformations 
and  refer  the  reader  to  the  references  cited  previously. 

2.5  Numerical  Procedure 

We  employ  a  numerical  procedure  similar  to  that  which  we  have  used  in  previous  work.  In 
particular,  a  rectangular  computational  grid  with  variable  grid  spacing  in  both  directions  is 
constructed.  The  spacing  of  grid  points  is  chosen  so  that  when  the  mapping  to  the  physical  space 
is  effected,  the  grid  will  be  small  in  regions  for  which  we  expect  a  strong  spatial  variation  and 
larger  in  regions  where  we  expect  a  more  gradual  variation.  The  derivatives  in  equations  22  and 
23  are  then  represented  by  standard  finite  differences,  and  the  equations,  which  may  be 
nonlinear,  are  solved  by  a  modified  Newton-Raphson  method  (Ames,  1977). 


3.  Calculations 


We  now  employ  the  previous  model  and  formalism  to  undertake  one  reasonably  simple 
calculation.  The  calculation  is  concerned  with  electrical  conduction  in  a  segment  of  the  rails  of  a 
railgun.  The  segment  is  assumed  to  be  sufficiently  far  behind  the  projectile  that  the  current  /  is 
conducted  only  in  the  z  direction.  If  we  view  a  cross  section  in  the  x-y  plane,  the  configuration 
appears  as  shown  in  figure  1 .  The  rail  on  the  left  conducts  current  in  the  positive  z  direction, 
while  the  rail  on  the  right  provides  the  return  current  path.  The  rails  have  height  h,  separation  w, 
and  thickness  s,  given  by  34,  44,  and  19.2  mm,  respectively.  It  is  evident  from  the  figure  that 
there  is  symmetry  about  y  =  0 ,  as  well  as  about  the  vertical  centerline  indicated.  Consequently, 

we  need  only  consider  one  quadrant  of  the  configuration  shown  and  choose,  for  convenience,  the 
upper  left-hand  quadrant. 

Shown  in  figure  2  are  the  grid  and  domain  D  for  the  calculation.  The  grid  shown  is  the  grid  in 
the  physical  space,  but  it  is  equivalent  to  the  computational  grid  because  of  the  transformation 
indicated  in  equation  24.  The  origin  of  the  x  axis  was  chosen  to  lie  at  the  left-hand  side  of  the 
domain.  Also  shown  in  the  figure  are  the  boundary  conditions  that  are  assumed  to  hold  on  the 
edges  of  the  domain.  The  top  and  left-hand  sides  are  assumed  to  be  sufficiently  far  from  the  rails 
that  we  can  assume  that  A  -  0.  From  symmetry  considerations,  we  also  have  A  —  0  along  the 
right-hand  side  (the  centerline  in  figure  1)  and  dA  /  dy  =  0  along  y  =  0. 
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For  simplicity,  we  had  intended  to  impose  a  constant  total  current  of  7o  for  all  calculations. 
Flowever,  in  order  to  avoid  discontinuities  in  either  the  current  or  its  derivatives,  we  assumed  a 
variation  of  the  form 


/  =  IQ  tanh(t/t0)  (25) 

and  chose  to  be  small  when  compared  to  the  total  calculation  time.  The  asymptotic  value  of  this 
functional  form  is  h,  and  I  / 10  is  greater  than  0.95  when  t  / 10  =2. 

It  is  then  necessary  to  determine  the  relationship  between  the  total  current  /  and  the  parameter  Eo 
that  appears  in  equation  10.  The  most  straightforward  way  is  to  integrate  equation  13  over  the 
cross-sectional  area  of  the  rail  and  equate  the  result  to  the  total  current  I.  A  somewhat  simpler 
expression  can  be  obtained,  however,  if  we  note  from  equation  10  that  the  current  density  Jean 
also  be  written  as 
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J  =  — -V  »VA, 


(26) 


where  we  have  assumed  that  the  permeability  is  constant  in  the  conductors.  If  we  now  exploit 
the  divergence  theorem  to  evaluate  the  required  integral,  we  have 


where  dA  /  dn  denotes  the  normal  derivative  of  A  at  the  boundary,  and  where  the  integral  is  to  be 
evaluated  over  the  boundaries  of  the  rail  cross  section.  The  procedure  is  standard  and  has  been 
discussed,  for  example,  by  Hoffmann  (1989).  The  factor  of  2  appears  on  the  right-hand  side 
because  only  the  top  half  of  the  rail  is  included  in  the  domain.  By  using  the  divergence  theorem, 
we  have  converted  a  2-D  integral  to  a  series  of  three  one-dimensional  integrals.  Simpson’s  rule 
is  used  to  perform  the  integrations. 

The  constraint  imposed  by  equation  27  should  be  sufficient  to  determine  E0.  It  has  proved, 
however,  to  be  enonnously  difficult  to  implement  this  relationship  in  our  computer  code  and  to 
obtain  accurate,  stable  values  for  E0.  Ultimately,  we  have  resorted  to  the  time-consuming 
procedure  of  assuming  a  value  for  E0  at  some  particular  time  step,  obtaining  convergence  of 
equations  10  and  14,  and  then  calculating  a  value  for  the  total  current,  denoted  by  Ic.  We  then 
changed  E0  and  recomputed.  This  process  allowed  us  to  find  the  acceptable  value  of  E0 — i.e., 
the  one  that  produced/  =  Ic — by  the  procedure  of  bracketing  and  bisection  (Press  et  ah,  1986). 
Even  then,  however,  we  sometimes  found  small  oscillations  in  Eq  as  a  function  of  time.  When 
these  small  oscillations  occurred,  we  fit  E0  to  a  curve,  typically  represented  by  a  three-parameter 
exponential  of  the  form 

E0  =  a  +  be~ct ,  (28) 

where  a,  b,  and  c  are  constants,  and  recomputed.  Usually,  we  could  then  obtain  reasonable 
agreement  between  the  computed  and  imposed  current.  This  laborious  procedure  is  not  entirely 
satisfactory,  but  we  have  been  unable  thus  far  to  produce  something  more  suitable.  The 
difficulty  of  obtaining  E0  for  situations  in  which  a  total  current  is  to  be  specified  has  been 
experienced  and  discussed  by  others  (Freeman,  1988). 

We  have  performed  calculations  for  a  number  of  different  values  of  h  as  well  as  for  a  number  of 
different  conditions.  We  discuss  only  one,  however.  For  this  particular  calculation  we  chose 
I0  =  500  kA,  t0  =100  |is,  and  assumed  that  the  imposed  current  was  represented  by  equation  25. 

The  temperature-dependent  thermophysical  properties  for  the  copper  rails,  shown  in  table  1, 
were  employed.  The  permeability  //  was  assumed  to  be  constant  and  equal  to  its  free-space  value 
/i o  everywhere. 
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We  followed  the  previously  discussed  procedure  to  determine  the  value  of  E 0  by  fitting  the 
numerical  results  to  an  exponential.  Approximate  values  of  the  constants  were  given  by: 
a  =  -16.3  V/m,  b  =  -2.14xl03  V/m,  and  c  =  1.34  x  104  s'1. 


Table  1 .  Properties  for  copper  employed  in  calculations 
(SI  units  are  used  unless  otherwise  noted). 


Property 

Value 

C 

360.0  +  o.  ir 

K 

412.0-0.078877 

p 

8900.0 

a"1  (gQ-cm) 

-0.542  +  0.00781T 

The  calculation  was  then  repeated  with  this  value  of  E o  and  produced  the  results  for  the 
computed  current  Ic  shown  in  figure  3;  the  actual  current  I,  given  by  equation  25,  is  also  shown 
for  comparison.  The  agreement  is  reasonable. 


Figure  3.  Actual  and  calculated  current  as  a  function  of  time. 
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The  first  results  we  show  are  the  vector  potential  contours  at  the  termination  of  the  calculation  at 
5  ms.  These  contours  are  shown  for  the  domain  D  in  figure  4;  they  can  be  sequentially  reflected 
about  the  lines  indicated  to  produce  somewhat  obvious  results  in  all  four  quadrants.  It  is  possible 
to  prove  (Powell,  2007)  that  the  vector  B  is  tangent  to  lines  of  constant  A.  Thus  the  arrows  in  the 
figure  indicate  the  direction  of  the  magnetic  induction  at  the  points  indicated.  The  magnitude  of 
By  varies  inversely  with  the  separation  of  the  contour  lines  in  the  x  direction,  whereas  the 
magnitude  of  Bx  varies  inversely  with  the  separation  in  the  y  direction.  Consequently,  the  B  field 
“wraps  around”  the  rails,  but  it  is  larger  within  the  bore  than  it  is  outside.  This  effect  is  a 
consequence  of  finite  rail  height  and  was  not  seen  in  our  previous  models  for  which  the  rail 
height  was  infinite.  It  is  also  of  some  interest  to  compare  the  maximum  value  of  Bv  in  the  bore 
with  the  value  obtained  for  an  infinite  rail  height  model.  The  maximum  value  of  Bv  occurs  at 
y  =  0  and  at  the  inner  surface  of  the  rail  (x  ~  22  cm);  its  value  there  is  about  8.2  T.  For  the 
infinite  rail  height  case,  Br  is  constant  in  the  bore  and  its  value  is  given  by  Bv  =  jij,  where  j 

represents  the  current  per  unit  rail  height.  For  j  ~  15  MA/m,  we  find  Br  ~  19  T  for  the  infinite 
rail  height  case.  There  is,  therefore,  a  difference  of  about  a  factor  of  2  for  the  two  cases. 


We  next  show  results  for  the  temperature  and  current  density  within  the  rail  at  several  times 
during  the  calculation.  The  current  density  is  shown  in  figure  5.  As  expected,  current  begins  to 
flow  just  inside  the  boundaries  of  the  conductor  where  the  gradients  in  the  magnetic  induction 
are  largest.  The  magnitude  is  particularly  high  near  the  corners  where  both  dB v  /  dy  and 
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dB v  /  dx  are  large  initially.  On  average,  the  magnitude  is  also  higher  on  the  right-hand  side  of  the 
rail  than  it  is  on  the  left-hand  side.  This  difference  results  because  of  the  proximity  of  the  return 
current  rail  on  the  right.  As  time  progresses,  the  current  diffuses  inward  from  the  sides  and  from 
the  top  of  the  rail.  At  late  times,  the  current  becomes  fairly  unifonn  within  the  rail  as  diffusion 
nears  completion.  The  magnitude  of  the  current  density  at  these  late  times  is  near  80  kA/cm2, 
i.e.,  the  total  current  divided  by  the  cross-sectional  area  of  the  rail.  An  approximate  time  for 
diffusion  to  occur  can  be  obtained  from  the  standard  skin-depth  formula,  S0  =  (mK/lkj))'  1 .  If 
we  assume  that  So  =  1  cm,  i.e.,  the  half-thickness  of  the  rail,  and  take  cr'1  to  be  its  room- 
temperature  value  of  about  1.8  pQ-cm,  we  obtain  a  time  of  about  2  ms.  That  result  seems 
consistent  with  results  in  figure  5. 
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Figure  5.  Current  density  contours  at  different  times. 
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The  isotherms,  shown  in  figure  6,  look  somewhat  similar  to  the  current  density  contours.  As  the 
current  diffuses  inward,  the  rails  are  heated  ohmically,  but  temperature  increases  lag  the  increase 
in  current  density.  Furthermore,  the  temperature  at  any  given  point  appears  to  rise  monotonically 
with  time.  This  behavior  indicates  that  ohmic  heating  is  much  more  significant  than  thermal 
conduction,  an  effect  that  is  small  both  because  of  the  relatively  short  time  scale  and  the  small 
temperature  gradients.  It  is  perhaps  surprising  that  the  temperature  rise  at  the  surface  of  the  rails 
(only  a  few  tens  of  degrees)  is  so  small.  Our  previous  calculations  with  the  infinite  rail  height 
models  have  indicated  substantially  higher  temperatures  in  the  rails.  The  difference  here  can  be 
attributed  to  two  factors.  First,  the  current  per  unit  rail  height  in  this  case  is  smaller  than  in  most 
of  our  previous  work,  and  the  initial  temperature  rise  at  the  rail  surface  varies  roughly  as  j  . 
Second,  some  significant  rail  heating  occurs  because  of  the  velocity  skin  effect  during  the  time 
that  the  trailing  edge  of  the  armature  is  near  a  given  point  on  the  rail  surface.  That  effect  is 
accounted  for  in  our  previous  calculations  but  not  in  this  one. 
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Figure  6.  Temperature  contours  at  different  times. 


We  have  finally  undertaken  an  approximate  calculation  of  the  inductance  gradient  L  for  this 
geometry  at  several  times.  This  quantity  can  be  represented  as  (Knoepfel,  1970) 
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(29) 


L'  =  -4r  f  f  B2  (x,  y)dxdy, 

fli 

where  the  integral  is  to  be  extended  over  the  entire  domain  D.  Values  for  L'  at  various  times 
(times  correspond  to  those  in  figures  3  and  4)  are  shown  in  table  2.  As  might  be  expected,  L' 
increases  as  the  current  diffuses  inward  from  the  sides  of  the  rail.  Kerrisk  (1981)  has  presented 
formulas  for  L'  in  the  high-frequency  limit  or,  equivalently,  at  very  early  times  when  the  current 
is  still  confined  to  the  rail  surfaces.  For  the  present  configuration,  Kerrisk’s  formulas  predict 
L'  =  0.56  |iH/m.  Similarly,  Grover  (1946)  has  obtained  formulas  that  apply  for  late-time,  steady- 
state  conditions  when  the  current  is  uniform  throughout  the  rail  cross  section.  His  results  predict 
L'  =  0.67  pH/m.  Results  in  table  2  lie  between  these  two  extreme  values,  as  should  be  expected. 
We  should  point  out  that  we  did  not  obtain  results  that  approached  Kerrisk’s  value  at  very  early 
times.  We  believe,  however,  that  this  discrepancy  arises  from  inadequate  resolution  of  the 
current  density  near  the  corners  when  the  current  is  largely  confined  to  that  location.  Indeed,  if 
we  decreased  the  grid  spacing  in  these  regions,  more  satisfactory  results  were  obtained.  This 
problem  will  be  addressed  in  more  detail  in  subsequent  work. 

Table  2.  Inductance  per  unit  length 
at  various  times. 


rfirs) 

L'  (uH/m) 

200 

0.591 

500 

0.598 

1000 

0.615 

1500 

0.627 

5000 

0.650 

4.  Summary  and  Future  Work 


We  have  extended  a  previously  developed  2-D  model  to  account  for  the  situation  in  which  there 
are  two  components  of  the  magnetic  induction  and  one  component  of  the  current  density.  The 
basic  formalism  was  established,  and  the  model  was  used  to  undertake  one  simple  calculation. 
The  calculation  is  concerned  with  the  diffusion  of  current  into  a  pair  of  infinitely  long  parallel 
rails.  The  configuration  approximates  diffusion  in  the  rails  of  a  railgun  at  distances  sufficiently 
far  behind  the  armature. 

In  future  efforts,  this  model  could  be  used,  for  example,  to  investigate  the  launch  or  deceleration 
of  metallic  plates  by  EM  fields  or  rail  erosion  in  railguns.  For  the  latter  problem,  it  will  likely  be 
necessary  to  include  the  trailing  edge  of  the  armature  and  to  account  in  some  approximate  way 
for  localized  heating  caused  by  the  velocity  skin  effect.  Such  a  calculation  will  extend  the  model 
beyond  its  rigorous  range  of  applicability,  but  it  may  offer  some  insight  into  how  erosion  occurs. 
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875  N  RANDOLPH  ST 
ARLINGTON  VA  22203-1995 

2  INST  FOR  ADVNCD  TCHNLGY 
THE  UNIV  OF  TEXAS 

AT  AUSTIN 
K  HSIEH 
I  MCNAB 

3925  W  BRAKER  LN 
AUSTIN  TX  78759-5316 

1  NORTH  CAROLINA 
STATE  UNIV 

DEPT  OF  NUCLEAR  ENGR 
M  BOURHAM 
BOX  7909 

RALEIGH  NC  27695-7909 

2  US  ARMY  TACOM  ARDEC 
FSAE  GCSS  TMA 

J  BENNETT 
D  LADD 
BLDG  354 

PICATINNY  ARSENAL  NJ  07806-5000 


ABERDEEN  PROVING  GROUND 


24  DIR  USARL 

AMSRD  ARL  WM  BA 
G  KATULKA 
AMSRD  ARL  WM  BD 
M  DEL  GUERCIO 
B  FORCH 
A  MICHLIN 
A  ZIELINSKI  (5  CPS) 
AMSRD  ARL  WM  MB 
J  TZENG 

AMSRD  ARL  WM  TA 
M  KEELE 

AMSRD  ARL  WM  TE 
P  BERNING 
C  HUMMER 
T  KOTTKE 
J  POWELL  (5  CPS) 

B  RINGERS 
E  SCHMIDT 
G  THOMSON 
C  UHLIG 
G  VUNNI 
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Intentionally  left  blank. 
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